Wolkewitz et al. Critical Care 201 4, 1 8:R64 
http://ccforum.eom/content/18/2/R64 



CRITICAL CARE 



RESEARCH Open Access 



Multilevel competing risk models to evaluate 
the risk of nosocomial infection 

Martin Wolkewitz 1,2 *, Ben S Cooper 3 , Mercedes Palomar-Martinez 4 , Francisco Alvarez-Lerma 5 , 
Pedro Olaechea-Astigarraga 6 , Adrian G Barnett 7 , Stephan Harbarth 8 and Martin Schumacher 1 



Abstract 

Introduction: Risk factor analyses for nosocomial infections (NIs) are complex. First, due to competing events for Nl, 
the association between risk factors of Nl as measured using hazard rates may not coincide with the association using 
cumulative probability (risk). Second, patients from the same intensive care unit (ICU) who share the same 
environmental exposure are likely to be more similar with regard to risk factors predisposing to a Nl than patients from 
different ICUs. We aimed to develop an analytical approach to account for both features and to use it to evaluate 
associations between patient- and ICU-level characteristics with both rates of Nl and competing risks and with the 
cumulative probability of infection. 

Methods: We considered a multicenter database of 1 59 intensive care units containing 1 09,21 6 admissions (81 3,739 
admission-days) from the Spanish HELICS-ENVIN ICU network. We analyzed the data using two models: an etiologic 
model (rate based) and a predictive model (risk based). In both models, random effects (shared frailties) were 
introduced to assess heterogeneity. Death and discharge without Nl are treated as competing events for Nl. 
Results: There was a large heterogeneity across ICUs in Nl hazard rates, which remained after accounting for 
multilevel risk factors, meaning that there are remaining unobserved ICU-specific factors that influence Nl occurrence. 
Heterogeneity across ICUs in terms of cumulative probability of Nl was even more pronounced. Several risk factors 
had markedly different associations in the rate-based and risk-based models. For some, the associations differed in 
magnitude. For example, high Acute Physiology and Chronic Health Evaluation II (APACHE II) scores were associated 
with modest increases in the rate of nosocomial bacteremia, but large increases in the risk. Others differed in sign, for 
example respiratory vs cardiovascular diagnostic categories were associated with a reduced rate of nosocomial 
bacteremia, but an increased risk. 

Conclusions: A combination of competing risks and multilevel models is required to understand direct and indirect 
risk factors for Nl and distinguish patient-level from ICU-level factors. 



Introduction 

Nosocomial infections (NIs) are a major threat for hospi- 
talized patients, particularly in intensive care units (ICUs), 
because they are associated with increased mortality and 
morbidity [1,2]. Analysis of data from large multicenter 
studies has the potential to improve our understanding 
of how patient- and ICU-level characteristics impact Nl 
outcomes. Such analysis is, however, complicated by two 
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factors: unexplained ICU-level variation and the impor- 
tance of competing risks. 

First, there are endogenous and exogenous modes of 
Nl acquisition [3]. Harbarth et al. [3] showed that about 
20% of NIs are exogenous and therefore potentially pre- 
ventable. Potential transmission routes of exogenous NIs 
are contact with contaminated environmental surfaces or 
cross-transmission via health-care workers or patients. 
Thus, patients from the same ICU who share the same 
environmental exposure are likely to be more similar 
with regard to acquiring a Nl than patients from differ- 
ent ICUs. In addition to patient-individual characteristics, 
ICU-specific factors, such as number of beds and nurses, 
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type of ICU and infection control policies, are poten- 
tial determinants for the occurrence of NI. To distinguish 
patient-level and ICU-specific factors requires multilevel 
analysis, but this is rarely used in hospital epidemiology 
[4]. The clustered structure of the data (patients within 
ICUs) often contains information that can be of value in 
understanding associations between risks and NIs [5]. 

Second, the length of ICU stay is a key determinant of 
the risk of NI. However, most patients are discharged from 
ICU or die in ICU without NI. Factors that are associ- 
ated with a high increased rate of infection are often also 
associated with an increased risk of dying in the hospital 
as well as with an extended length of stay. These compet- 
ing events play an important role in risk interpretation of 
NI and make extended survival models necessary [6-9]. 
Again, rates of NI, discharge and death without NI might 
also depend on patient- as well as on ICU-level factors. 
Thus, a combination of extended survival and multilevel 
models is required to understand how different risk fac- 
tors impact NI outcomes. The aim of this paper is to apply 
established and innovative statistical methods [10,11] to 
investigate heterogeneity in risks and rates across ICUs of 
NI and concurrently occurring competing events. 

Material and methods 

Spanish ICU data 

We used a multicenter database from the Spanish surveil- 
lance network HELICS-ENVIN [12], embedded in the 
HELICS project (Hospitals in Europe Link for Infection 
Control through Surveillance). The reliability and quality 
of the surveillance program has recently been investigated 
[13]. Data were prospectively collected on an individ- 
ual patient level and also aggregate ICU level. For our 
purpose, we included ICUs that contributed to the reg- 
istry between January 2006 and December 2011 and we 
included only patients who stayed at least two days in ICU. 
We excluded ICUs that contributed less than 100 patient 
admissions to the cohort to reduce artificial heterogene- 
ity. To get a robust outcome, we focused on primary or 
secondary nosocomial bacteremia (NB). The study pop- 
ulation, 159 intensive care units with 109,216 admissions 
(813,739 admission-days), is summarized in Table 1. The 
data of this official surveillance are encrypted and com- 
pletely anonymous. Patients' consent was not needed. 
This Deutsche Forschungsgemeinschaft research project 
was approved by the ethics committee of the University 
Medical Center Freiburg, Germany. 

Shared frailty models for competing risks 
Model 1: etiologic model (rate-based) 

The classic way to analyze competing risks data is to study 
event-specific hazard rates, i.e., fitting a proportional 
hazard model for each event (NB, death without NB and 
discharge without NB) separately. Random effects for each 



ICU (i.e. frailties) can be introduced by a shared gamma 
frailty model [10,14] (see details in Additional file 1). For 
each of the three events (NB, death and discharge), we fit- 
ted models to assess heterogeneity for NB rates, quantified 
by the corresponding variance estimator 9. Large vari- 
ances signify a closer similarity between patients within 
ICU and greater heterogeneity across ICUs. The follow- 
ing quantities were calculated for each of the three events: 
the baseline hazard, ICU effects, variance of ICU effects 
(9) and the hazard ratios for multilevel risk factors at 
the patient and ICU levels. In this approach, the hazards, 
i.e., the daily risks of the primary outcome (NB) and the 
competing events (death or discharge without NB), are 
studied. Note that the hazard of NB does not depend on 
the competing events. 

Model 2: predictive model (risk-based) 

The cumulative incidence function of NB is defined as the 
probability of NB over a period of time and interpreted 
as the actual risk of NB occurring in this time period. 
This approach is useful for predicting NB. It has previ- 
ously been shown that the way in which risk factors are 
associated with the NB hazard (instantaneous risk) may 
not coincide with the way these factors are associated with 
the cumulative incidence of NB (cumulative risk) [15]. 
To study the risk (cumulative incidence function) of NB, 
we used the Fine and Gray model [16] and introduced a 
shared frailty structure to investigate heterogeneity in a 
similar way as Katsahian et al. [11]. Using this model we 
calculated the following quantities: the baseline subdistri- 
bution hazard and corresponding cumulative incidence of 
NB, ICU effects, variance of ICU effects (9) and the sub- 
distribution hazard ratios for multilevel risk factors at the 
patient and ICU levels. In this approach, the cumulative 
risk of the primary outcome (NB) is studied. In contrast 
to the event-specific approach, the cumulative risk of NB 
depends on the NB hazard as well as on the hazards of the 
competing events and it tends towards the overall risk, i.e., 
the incidence proportion of NB. 

For both models, we first used a model with frail- 
ties for each ICU but without covariates (null model). 
Then, we considered a multivariate model by introducing 
patient-individual as well as ICU-specific covariates and 
estimated the frailties for each ICU. For all analyses we 
used the flexible R package frailtyPack [17]. 

Results 

In the following we present the detailed results for pri- 
mary and secondary nosocomial bacteremia (NB). 

Baseline hazard rates and cumulative incidence function 

The overall baseline hazard rates based on the null models 
without covariates are shown for each event in Figure 1. 
The hazard rate of NB is increasing with the time from 
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Table 1 Description of study population 

General 

Number of admissions 
Number of admission-days 
Number of ICUs 

Number of nosocomial bacteremia during ICU stay 
Number of deaths without NB during ICU stay 
Number of discharges without NBfrom ICU 
Number of administrative censored admissions 
Overall risk of nosocomial bacteremia (censored excluded) 
Overall rate of nosocomial bacteremia 



Frequency 

109,216 
813,739 
159 

5,498 (5.03%) 
12,678(11.61%) 
90,142 (82.54%) 
898 (0.82%) 
5.08% 

6.75 / 1 ,000 admission-days 



Risk factors 

ICU / hospital level covariates 

Number of beds in ICU: 

0 to 1 0 (reference) 

1 1 to 20 

21 to 30 

31 to 40 

>40 

Number of beds in hospital: 
0 to 500 (reference) 
501 to 1,000 
> 1,000 

Type of hospital: 

Private 

Public 

Type of ICU: 

Polyvalent (reference) 

Medical 

Surgery 

Coronary 

Traumatology 

Post-surgery cardiology 

Burn 

University+teaching hospital (reference) 
Teaching hospital (no university) 
Hospital without teaching/university 

Calendar year of admission 

2006 (reference) 

2007 

2008 

2009 

2010+ 



Frequency (%) 
(patient level) 



30,389 (27.82) 
46,524 (42.60) 
21,668(19.84) 
7,673 (7.03) 
2,962 (2.71) 

52,426 (48.00) 
8,259 (7.56) 
48,531 (44.44) 

6,541 (5.99) 
102,675 (94.01) 

96,478 (88.34) 
3,895 (3.57) 
3,103 (2.84) 
525 (0.48) 
3,273 (3.00) 

I, 834(1.68) 
108 (0.10) 

67,917(62.19) 
30,089 (27.55) 

II, 210(10.26) 

14,318(13.11) 
17,819(16.32) 
21,559(19.74) 
25,660 (23.49) 
29,860 (27.34) 



Frequency (%) 
(ICU level) 



61 (38.36) 
71 (44.65) 
19(11.95) 
5(3.14) 
3 (1.89) 

94 (59.12) 
1 0 (6.29) 
55 (34.59) 

12(7.55) 
147 (92.45) 

1 38 (86.79) 
5(3.14) 

3 (1.89) 
2(1.26) 
6 (3.77) 

4 (2.52) 
1 (0.63) 

86 (54.09) 
52 (32.70) 
21 (13.21) 
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Table 1 Description of study population (Continued) 



Patient level covariates 

APACHE II score: 
0 to 1 0 (reference) 
1 1 to 20 
21 to 30 
>30 

Age (years): 
Oto 40 
41 to 60 

61 to 80 (reference) 



Days in hospital before ICU admission: 
0 to 3 (reference) 
4 to 6 
7 to 10 
> 10 

Type of diagnosis: 
Cardiovascular (reference) 
Respiratory 
Gastrointestinal 
Central nervous system 
Other diagnoses 

Antibiotic treatment 48 h before and/or after ICU admission 
Gender (male) 

Origin: community (reference) 
Origin: hospital/ICU 
Trauma 



40,353 (36.95) 
44,654 (40.89) 
19,191 (17.57) 
5,018(4.59) 

54,477 (49.88) 
12,931 (11.84) 
29,371 (26.89) 
12,437 (11.39) 

87,208 (79.85) 
5,864 (5.37) 
4,962 (4.54) 
11,182(10.24) 



54,374 (49.79) 
15,243 (13.96) 
14,626(13.39) 
17,567 (16.08) 
7,406 (6.78) 

23,178 (21.22) 
71,223 (65.21) 
54,996 (50.36) 
54,220 (49.64) 
8,927(8.17) 
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Figure 2 Estimated subdistributor! hazard function (left) and cumulative incidence function (right). Data are from the shared frailty model 
for the subdistribution hazard of NB without covariates. The subdistribution hazard is shown as the black curves and associated 95% confidence 
intervals as broken lines. 



admission and has a peak at day 15; for instance, the daily 
risk of acquiring a NB at day 10 (day 30) for a patient is 
about 1% (1.5%) given that he or she has stayed at ICU 
without a NB for at least 9 days (29 days). The death haz- 
ard rate without NB is about 2% for the first 40 days from 
admission. Obviously, the discharge hazard rate without 
NB is the strongest hazard with its peak (about 25%) about 
5 to 6 days after admission and a strong decrease after- 
wards; meaning the likelihood of discharge without NB 
decreases for each survived day in ICU. The subdistribu- 
tion hazard function and cumulative incidence function 
of NB are displayed in Figure 2; the cumulative incidence 
function tends towards the overall risk of 5.08%. 

Heterogeneity across ICUs 

The corresponding ICU effects (random effects or frail- 
ties) of the null models are shown in Figure 3; for instance, 
an ICU effect of 2 means that the baseline hazard of this 
ICU is twice as large as the hazard averaged over all ICUs. 



The observed heterogeneity across ICUs in the rates of 
NB is remarkably large (9 = 0.26 with standard error 
(SE) = 0.038), in contrast to the heterogeneity in rates of 
death without NB (9 = 0.14 with SE = 0.019) and discharge 
without NB (9 = 0.15 with SE = 0.017). The heterogeneity 
in risks of NB is even larger than in the rates (estimated by 
9 = 0.64 with SE = 0.076); see Figure 4. Large heterogene- 
ity means that event times are strongly correlated within 
ICUs. One might conclude that observed and unob- 
served ICU-specific factors play a more substantial role 
in the NB hazard rates compared to the competing event 
rates. 

Including covariates in the models led to only a small 
reduction in ICU-level heterogeneities for all three of the 
competing outcomes, with 9 falling to 0.19 for NB, still 
larger than the heterogeneities in rates of the competing 
events (Table 2). Separate models (only ICU-level factors 
and only patient-level factors) showed that ICU-level fac- 
tors reduced heterogeneity more than the patient-level 
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Figure 3 Estimated frailties (random effects) for each ICU for the three outcomes. Data are from the null model without covariates. The ICUs 
are ordered on the X-axis according to the contributing number of patients; the circles are proportional to the magnitude of contribution of the ICU, 
e.g., an ICU with sqrt(Npatients) = 60 contributed 60 2 = 3600 admissions to the cohort. The V-axis has a log scale. NB, nosocomial bacteremia. 
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Figure 4 Estimated frailties (random effects) for each ICU. Data are from the shared frailty model for the subdistribution hazard of NB without 
covariates. The ICUs are ordered on the X-axis according to the contributing number of patients; the circles are proportional to the magnitude of 
contribution of the ICU, e.g., an ICU with sqrt(Npatients) = 60 contributed 60 2 = 3600 admissions to the cohort. The /-axis has a log scale. 



factors (Table 2). It follows that the impact of unobserved 
ICU-specific factors is large. The heterogeneity reduction 
was stronger in the NB risk than the NB rate model. There 
are patient-level factors, such as the APACHE II score, 
type of diagnosis and trauma (Table 2), that have a strong 
effect in reducing the discharge without NB hazard (the 
strongest competing risk hazard), i.e., patients with high 
APACHE II scores or trauma stay longer in ICU. There are 
also ICUTevel factors that are more pronounced: num- 
ber of beds in ICU and teaching hospital. Thus, there is 
an indirect effect on the risk for NB leading to increased 
subdistribution hazard ratios. And strong effects lead to a 
reduction in heterogeneity. 

Multilevel risk factors: patient level 

Results from the risk factor analysis on the patient as 
well as ICU level are shown in Table 2. The hazard ratios 
for NB from the event-specific analysis (model 1) reflect 
direct effects on the occurrence of NB; this event-specific 
approach also shows how these factors are associated 
with the competing events, i.e., those which have poten- 
tially indirect effects on NB. The subdistribution analysis 
(model 2) is a summary analysis and studies the effects 
on the risk of NB (as a cumulative incidence function). 
Factors, i.e., an exposure, with a hazard ratio lower than 
1 for the competing events increase the risk of NB since 
exposed patients stay longer in the ICU leading to more 
NBs in the exposed group. For instance, traumatic patients 
are associated with an increased daily risk of NB (hazard 



rate = 1.13) and also with a decreased daily risk of being 
discharged or dying without NB (both hazard rates lower 
than 1). In other words: traumatic patients acquire NB 
more frequently per day and in addition, as they stay 
longer at ICU, their extended time at ICU is itself a risk 
factor for NB, which causes them to acquire more NBs 
eventually. This indirectly adds to the cumulative risk of 
NB (subdistribution hazard rate = 1.81). 

In a similar way, factors with a hazard ratio greater 
than 1 for the competing events decrease the risk of NB. 
For example, patients older than 80 years acquire fewer 
NBs per day than patients aged 61 to 80 years (hazard 
rate = 0.77) but they also die (and are discharged) 
more frequently per day at ICU (both hazard rates 
are greater than 1). This indirect effect additionally 
reduces the cumulative risk of NB (subdistribution hazard 
rate = 0.51). Often, the hazard ratios for discharge and 
death without NB are diametrically opposed; then, the 
discharge hazard has usually a stronger indirect effect 
because of its larger magnitude compared with the death 
hazard (Figure 1). For instance, the APACHE II score is 
highly associated with an increased rate of NB but also 
highly associated with an increased death (without NB) 
rate and a decreased discharge (without NB) rate (i.e., 
with a longer ICU stay). The indirect effect due to an 
extended length of ICU stay makes patients with higher 
APACHE II score acquire an NB more frequently, which 
is quantified by the subdistribution hazard ratios. It is 
even possible for a hazard ratio for NB to be lower than 



Table 2 Results from multivariate analysis 

Risk factors 



Event-specific analysis 



Subdistribution analysis 



a 
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ICU / hospital-level covariates 

Number of beds in ICU: 
11-20 vs 0-10 
21-30 vs 0-10 
31-40 vs 0-10 
>40 vs 0-10 

Number of beds in hospital: 
> 1 000 vs 0-500 
501-1000 vs 0-500 

Teaching hospital (only) vs university+teaching 
No teaching hospital vs university+teaching 
Type of hospital (private vs public) 

Type of ICU: 

Medical vs polyvalent 
Surgery vs polyvalent 
Coronary vs polyvalent 
Traumatology vs polyvalent 
Post-surgery cardiology vs polyvalent 
Burn vs polyvalent 

Calendar year of admission 

2007 vs 2006 

2008 vs 2006 

2009 vs 2006 
2010+vs 2006 



Nosocomial bacteremia (NB) 

HR(95%CI) 



1 .1 8 (0.96-1 .46) 
1.31 (0.97-1.77) 
1.56(0.98-2.50) 
1.09 (0.55-2.18) 

1.30(0.85-1.99) 
1.11 (0.91-1.36) 

1 .00 (0.82-1 .22) 
0.71 (0.51-0.97) 
0.98 (0.70-1.39) 

0.71 (0.44-1.14) 
1.07 (0.62-1.85) 
0.71 (0.28-1.82) 
1.20(0.79-1.82) 
1.47(0.85-2.54) 
0.48(0.16-1.42) 

1.02 (0.92-1.13) 
0.98 (0.89-1.08) 
1.09 (0.99-1.19) 
0.84 (0.76-0.93) 



Death without NB in ICU 

HR (95% CI) 



1.01 (0.87-1.17) 
0.86 (0.69-1.08) 
0.81 (0.57-1.14) 

1.00 (0.60-1.65) 

0.94 (0.70-1.27) 
1.04 (0.90-1.21) 

1.03 (0.89-1.20) 
1.08 (0.86-1.34) 
0.82 (0.64-1.05) 

1.07 (0.78-1.47) 
0.60 (0.40-0.91) 
1.27 (0.64-2.51) 
1.03 (0.76-1.40) 
0.71 (0.46-1.08) 
0.70 (0.31-1.62) 

1.01 (0.94-1.08) 
1.09(1.02-1.16) 

1.02 (0.96-1.09) 
1.08(1.01-1.15) 



Discharge without NB 

HR (95% CI) 



1.43 (1.23-1.67) 
1.15 (0.89-1.47) 
1.83 (1.22-2.74) 
1 .37 (0.76-2.48) 

1.30 (0.90-1.87) 
1.13 (0.95-1.34) 

1.87(1.61-2.18) 
2.42(1.92-3.03) 
0.80 (0.61-1.05) 

1.44 (0.99-2.10) 
1.20 (0.73-1.96) 
0.83 (0.44-1.57) 
1.28 (0.90-1.81) 
1.95 (1.27-3.01) 
0.70 (0.31-1.59) 

0.95 (0.92-0.97) 
0.93 (0.91-0.95) 
0.86 (0.83-0.88) 
0.94(0.91-0.96) 



Nosocomial bacteremia 

Subdistribution HR (95% CI) 



1.34 (0.99-1.81) 
2.17 (1.42-3.30) 
1.85 (0.95-3.62) 
1.37 (0.50-3.72) 

1.32 (0.71-2.47) 
1.16(0.87-1.53) 

0.94 (0.71-1.26) 
0.59(0.38-0.91) 
1.18(0.76-1.85) 

0.57 (0.30-1.08) 
1.13 (0.51-2.50) 
0.73 (0.24-2.28) 
1.15(0.64-2.05) 
1.39(0.63-3.05) 
1.20 (0.30-4.79) 

1.02 (0.93-1.13) 

1.03 (0.93-1.14) 
1.19(1.09-1.31) 
0.83 (0.75-0.91) 



Table 2 Results from multivariate analysis (Continued) 
Patient level covariates 

APACHE II score: 

1 1-20 vs 0-10 1.30(1.20-1.42) 

21-30vs0-10 1.38 (1.26-1.50) 

> 31 vs 0-10 1.54(1.37-1.73) 

Age (years): 

0-40 vs 61 -80 1.08(0.99-1.18) 

40-60 vs 61-80 1.09(1.02-1.16) 

> 80 vs 6 1 -80 0.77 (0.68-0.86) 
Days in hospital before ICU admission: 

4-6 vs 0-3 1.12 (1.00-1.25) 

6-10vs0-3 1.14(1.01-1.28) 

>10vs0-3 1.17 (1.08-1.27) 
Type of diagnosis: 

Respiratory vs cardiovascular 0.86 (0.80-0.94) 

Gastrointestinal vs cardiovascular 1.18(1 .09-1 .29) 

Central nervous system vs cardiovascular 0.91 (0.84-0.99) 

Other diagnoses vs cardiovascular 1 .20 (1 .08-1 .32) 



1.99(1.85-2.14) 
4.11 (3.83-4.41) 
6.65 (6.15-7.19) 

0.59 (0.55-0.64) 
0.79 (0.76-0.83) 
1.65 (1.57-1.73) 

1.09(1.01-1.17) 
1.18(1.10-1.27) 
1.22 (1.16-1.29) 

0.93 (0.88-0.98) 
0.96 (0.91-1.02) 
1.41 (1.34-1.49) 
0.78(0.71-0.85) 



0.53 (0.52-0.54) 
0.28 (0.27-0.28) 
0.17(0.17-0.18) 

1.05 (1.03-1.07) 
0.98 (0.97-1.00) 
1.18(1.15-1.21) 

0.86 (0.83-0.89) 
0.85 (0.82-0.87) 
0.78 (0.76-0.80) 

0.63 (0.62-0.65) 
0.78 (0.76-0.79) 
0.65 (0.64-0.66) 
0.76 (0.74-0.78) 



2.81 (2.59-3.06) 
4.54(4.15-4.96) 
5.78 (5.14-6.49) 

1.12(1.02-1.22) 
1.18(1.11-1.26) 
0.51 (0.45-0.57) 

1.27(1.13-1.42) 
1.33 (1.18-1.49) 
1.43 (1.32-1.55) 

1.43 (1.32-1.55) 
1.65 (1.51-1.80) 
1.38(1.28-1.50) 
1.96(1.78-2.16) 



Antibiotic treatment 48 h before and/or after ICU admission 0.83 (0.77-0.88) 

Gender 1.08(1.01-1.14) 

Origin (hospital/ICU vs community) 1 .05 (0.98-1 .1 2) 

Trauma 1.13(1.04-1.24) 



1.08(1.04-1.13) 
1.01 (0.97-1.05) 
0.98 (0.94-1 .02) 
0.67 (0.62-0.72) 



0.70 (0.69-0.71) 
0.99 (0.98-1.01) 
1.00 (0.99-1.02) 
0.64 (0.62-0.66) 



1.01 (0.95-1.09) 
1.12(1.06-1.18) 
1.09(1.02-1.17) 
1.81 (1.66-1.98) 



Variance of heterogeneity (without covariates) 0.26 (SE 0.038) 

Variance of heterogeneity (with ICU-level covariates only) 0.20 (SE 0.030) 

Variance of heterogeneity (with patient-level covariates only) 0.25 (SE 0.037) 

Variance of heterogeneity (with all covariates) 0.1 9 (SE 0.030) 



0.14 (SE 0.019) 
0.1 2 (SE 0.01 5) 
0.1 2 (SE 0.01 6) 
0.11 (SE0.015) 



0.1 5 (SE 0.01 7) 
0.09 (SE 0.011) 
0.1 7 (SE 0.025) 
0.1 7 (SE 0.024) 



0.64 (SE 0.076) 
0.47 (SE 0.061) 
0.56 (SE 0.069) 
0.40 (SE 0.052) 



CI, confidence interval; HR, hazard rate; NB, nosocomial bacteremia; SE, standard error. 
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1 but for the subdistribution hazard ratio to be larger 
than 1 (such as respiratory vs cardiovascular diagnosis). 
This is due to an indirect effect on the competing events: 
patients with a respiratory diagnosis stay much longer 
in ICU without NB than patients with a cardiovascular 
diagnosis, or in other words, their daily risk of being dis- 
charged or dying without NB is reduced (hazard ratios 
of both competing events are lower than 1). This addi- 
tional impact on the risk of NB is so strong that respiratory 
patients acquire more NB than cardiovascular patients 
even though their daily risk of NB is reduced. These 
results highlight just some the complexities of exam- 
ining risk factors of NI in the presence of competing 
risks. 

Multilevel risk factors: ICU level 

A larger number of beds in an ICU was somewhat indi- 
rectly associated with an increased risk of NB (Table 2) 
because the number of beds in ICUs were moderately 
associated with the competing events. The rate of NB 
was higher in university/teaching hospitals compared 
with those without teaching, an effect which has also 
been found elsewhere [18]. This effect is even more pro- 
nounced in the subdistribution model because patients in 
non-teaching hospitals are discharged quicker. 

Discussion 

In this paper, we used two multilevel competing risks 
models to evaluate risk factors of NI. To our knowl- 
edge, this is the first study to investigate the heterogeneity 
across ICUs in risks and rates of NI using a large multicen- 
ter cohort accounting for both ICU clustering effects and 
competing risks. Further, we showed that it is necessary 
to perform a multilevel competing risk analysis to under- 
stand fully the direct and indirect effects of risk factors on 
the occurrence of NI. 

Our findings have the following implications. First, the 
large heterogeneity indicates that the impact of unob- 
served ICU-specific factors on the risk and rate of NI 
is large, even after accounting for important patient- 
and observed ICU-level characteristics. Thus, surveillance 
networks are encouraged to collect further potential ICU- 
level risk factors in addition to patient-level data. This 
large heterogeneity might also explain why ICU-based 
infection control strategies might work in some ICUs 
but not in others. It emphasizes the need for multicen- 
ter intervention studies rather than single-center pilot 
studies [19,20]. Ignoring heterogeneity in the analysis of 
multicenter studies can lead to biased results and mis- 
leading conclusions. In our cohort, the risks of NI were 
more heterogeneous than rates of NI across ICUs. From 
the mathematical point of view, this is not necessarily 
the case since different correlations between competing 



hazard constellations could potentially result in a very 
similar risk of NI. 

Second, competing risks play an essential role in the 
understanding of NI occurrence and the analysis must 
account for this [7]. The distinction between indirect and 
direct effects is a key issue for understanding the associa- 
tions between risk factors and NI outcomes. For instance, 
a cohort study of African children reported no associa- 
tion between burns and pediatric hospital-acquired bac- 
teremia [21] even though the cumulative risk is about 
three times higher, due to competing events [7]. There- 
fore, we recommend that the results from both mod- 
els (the rate and the risk models) are reported [22] to 
make hidden and indirect effects transparent. This is very 
important for NI since risk factors for NI are often also 
associated with the competing risks for NI as well. This 
may result in more pronounced effects in the risk model 
(such as APACHE II score, type of diagnosis or trauma) or 
effects that are relevant only in the rate model but irrel- 
evant in the risk model (e.g., antibiotic treatment before 
ICU admission). It could even lead to apparent diamet- 
rically opposed results (e.g., both respiratory and central 
nervous system vs cardiovascular diagnoses). Third, our 
approach has the potential to identify ICUs with unusu- 
ally high or low rates. Nevertheless, we emphasize that 
extreme caution is needed when using NI rates or risks for 
benchmarking or as quality indicators. Besides problems 
in reliability, validity [23] and statistical uncertainty, there 
are further complications from competing risks as shown 
here. 

There are some limitations to our study. As in other 
large surveillance studies based on volunteer ICUs, data 
can be subject to reporting, information or selection bias. 
Therefore, it is possible that part of the heterogeneity 
might be attributed to over- or under-reporting of NI 
rather than to real ICU factors [24]. Information bias 
might occur due to unreported NB cases. As suggested by 
Hansen et al. [18], surveillance data need to be validated 
to counter differences across ICUs. Further, an extrapola- 
tion of our findings to other European countries is limited 
due to differences in ICU management. 

Multilevel analyses in hospital infection epidemiology 
are still rare but necessary to evaluate effects of individual- 
level and group-level factors [4]. Vellinga et al. have pro- 
duced a valuable overview of the principles of multilevel 
analysis for antimicrobial resistance studies [5]. We have 
extended their methodology to complex survival data for 
hospital infections because modeling the timing of events 
(infection, death or discharge) is crucial. 

Conclusions 

We encourage further investigations using our method- 
ological approaches, e.g., to evaluate the occurrence of 
antimicrobial resistance by exploring antibiotic usage at 
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the patient level and ICU level simultaneously. Statis- 
tical models and corresponding codes are available in 
Additional file 1. 

Key messages 

• Discharge from and death in an intensive care unit 
are competing risks for nosocomial infection. 

• There are factors on the patient as well as on the 
ICU level influencing the occurrence of nosocomial 
infections. 

• Analysis of data from large multicenter studies has 
the potential to improve our understanding of how 
patient and ICU-level characteristics impact 
nosocomial infections. 

• A combination of multilevel and competing risk 
models are necessary to analyze such complex data. 

• We encourage further investigations by using our 
methodological approaches to evaluate our findings 
of unexplained heterogeneity. The statistical code is 
available in Additional file 1. 

Additional file 



Additional file 1 : Appendix: Statistical methods and R code for 
multilevel competing risk models. 
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